home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / JACOBI.CPP < prev    next >
C/C++ Source or Header  |  1995-01-11  |  3KB  |  100 lines
  1. //$$jacobi.cpp                           jacobi eigenvalue analysis
  2.  
  3. // Copyright (C) 1991,2,3,4: R B Davies
  4.  
  5.  
  6. //#define WANT_STREAM
  7.  
  8.  
  9. #define WANT_MATH
  10.  
  11. #include "include.h"
  12. #include "newmat.h"
  13. #include "precisio.h"
  14. #include "newmatrm.h"
  15.  
  16.  
  17.  
  18. void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A,
  19.    Matrix& V, Boolean eivec)
  20. {
  21.    Real epsilon = FloatingPointPrecision::Epsilon();
  22.    Tracer et("Jacobi");
  23.    int n = X.Nrows(); DiagonalMatrix B(n), Z(n); D.ReDimension(n); A = X;
  24.    if (eivec) { V.ReDimension(n,n); D = 1.0; V = D; }
  25.    B << A; D = B; Z = 0.0; A.Inject(Z);
  26.    for (int i=1; i<=50; i++)
  27.    {
  28.       Real sm=0.0; Real* a = A.Store(); int p = A.Storage();
  29.       while (p--) sm += fabs(*a++);            // have previously zeroed diags
  30.       if (sm==0.0) return;
  31.       Real tresh = (i<4) ? 0.2 * sm / square(n) : 0.0; a = A.Store();
  32.       for (p = 0; p < n; p++)
  33.       {
  34.      Real* ap1 = a + (p*(p+1))/2;
  35.          Real& zp = Z.element(p); Real& dp = D.element(p);
  36.      for (int q = p+1; q < n; q++)
  37.      {
  38.         Real* ap = ap1; Real* aq = a + (q*(q+1))/2;
  39.             Real& zq = Z.element(q); Real& dq = D.element(q);
  40.             Real& apq = A.element(q,p);
  41.             Real g = 100 * fabs(apq); Real adp = fabs(dp); Real adq = fabs(dq);
  42.  
  43.         if (i>4 && g < epsilon*adp && g < epsilon*adq) apq = 0.0;
  44.         else if (fabs(apq) > tresh)
  45.         {
  46.            Real t; Real h = dq - dp; Real ah = fabs(h);
  47.            if (g < epsilon*ah) t = apq / h;
  48.            else
  49.            {
  50.           Real theta = 0.5 * h / apq;
  51.           t = 1.0 / ( fabs(theta) + sqrt(1.0 + square(theta)) );
  52.           if (theta<0.0) t = -t;
  53.            }
  54.            Real c = 1.0 / sqrt(1.0 + square(t)); Real s = t * c;
  55.            Real tau = s / (1.0 + c); h = t * apq;
  56.                zp -= h; zq += h; dp -= h; dq += h; apq = 0.0;
  57.            int j = p;
  58.            while (j--)
  59.            {
  60.           g = *ap; h = *aq;
  61.           *ap++ = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
  62.            }
  63.            int ip = p+1; j = q-ip; ap += ip++; aq++;
  64.            while (j--)
  65.            {
  66.           g = *ap; h = *aq;
  67.                   *ap = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
  68.           ap += ip++;
  69.            }
  70.            int iq = q+1; j = n-iq; ap += ip++; aq += iq++;
  71.            while (j--)
  72.            {
  73.           g = *ap; h = *aq;
  74.           *ap = g-s*(h+g*tau); *aq = h+s*(g-h*tau);
  75.           ap += ip++; aq += iq++;
  76.            }
  77.            if (eivec)
  78.            {
  79.                   RectMatrixCol VP(V,p); RectMatrixCol VQ(V,q);
  80.                   Rotate(VP, VQ, tau, s);
  81.            }
  82.         }
  83.      }
  84.       }
  85.       B = B + Z; D = B; Z = 0.0;
  86.    }
  87.    Throw(ConvergenceException(X));
  88. }
  89.  
  90. void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D)
  91. { SymmetricMatrix A; Matrix V; Jacobi(X,D,A,V,FALSE); }
  92.  
  93. void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A)
  94. { Matrix V; Jacobi(X,D,A,V,FALSE); }
  95.  
  96. void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, Matrix& V)
  97. { SymmetricMatrix A; Jacobi(X,D,A,V,TRUE); }
  98.  
  99.  
  100.